!
! 	routine de la résolution de Au = g par la méthode d Euler
!
SUBROUTINE EULER(n,Nf,nu,h,dt,a,b,c,seuilResidu,nbIterMax,u)!
! 	Arguments
!   	 u : OUT : Vecteur résultat du système linéaire
!   	 n : IN : dimension des blocs dans la matrice A
!		 a, nu, b, h, c : IN : parametres qui determinent quelle equation est utilisee
!
implicit none
INTEGER n, Nf, k, T, I, J
! Vecteur u(k+1)
DOUBLE PRECISION, DIMENSION(n*n):: u
! Vecteur u(k)
DOUBLE PRECISION, DIMENSION(n*n):: uk
! Matrice A et g
DOUBLE PRECISION Mat(5), f(n*n)
DOUBLE PRECISION nu, a, b, c, h, dt
DOUBLE PRECISION D(1)
DOUBLE PRECISION RES(n*n), vect_residu(n*n)
DOUBLE PRECISION residu
DOUBLE PRECISION ti, xi, yj
DOUBLE PRECISION :: t_debut, t_fin
INTEGER compteur
! Critere d'arret
DOUBLE PRECISION, INTENT(IN) :: seuilResidu
INTEGER, INTENT(IN) :: nbIterMax
REAL :: pi
! On lance le chronometre
CALL CPU_TIME(t_debut)
! Calcul de Pi
pi = ACOS(0.0)
! Initialisation a 0
uk = 0*uk
residu = 1.0D0
!
! Appel des differentes subroutines
CALL createMatTemporelle(nu,a,b,h,Mat)
!
! Initialisation du compteur
compteur = 1
! On applique la methode d'euler
!
! A LAISSER, SINON CA MARCHE PLUS
PRINT *, ''
DO WHILE ((residu .GT. seuilResidu) .AND. (compteur .LT. nbIterMax))
	! Calcul pour chaque iteration de temps
	DO T = 1, 2
		ti = 1.0*T*dt
		k = 1
		DO I = 1, n
			! Discretisation du domaine
			! On force la division a etre reelle
			xi = 1.0*I*h 
			DO J = 1, n
				! On force la division a etre reelle egalement
				yj = 1.0*J*h
				! Creation de g
				! ATTENTION, c’est un vecteur colonne
				f(k) = pi*EXP(-pi*ti)*(-SIN(pi*xi)*SIN(pi*yj) + a*COS(pi*xi)*SIN(pi*yj) &
				& + b*SIN(pi*xi)*COS(pi*yj) - 2*nu*pi*SIN(pi*xi)*SIN(pi*yj))
				k = k+1
			 ENDDO
		ENDDO
		! On insere le calcul par la methode euler dans la construction de f
		! Calcul du produit matrice vecteur A*u
		CALL produitmatricevecteur(Mat,uk,n,Nf,RES)
		! Calcul de u(k+1) par la méthode de Jacobi
		! On utilise la subroutine combinaison lineaire, car sinon on obtient seulement un scalaire
		CALL COMBINAISONLINEAIRE(n*n,dt,f - RES,1.0,uk,u)
		WRITE(*,*) 'u',u
		! Calcul du residu
		CALL produitmatricevecteur(Mat,u,n,n,vect_residu)
		vect_residu = f - vect_residu
		! Calcul de la norme via le produit scalaire
		CALL PRODUITSCALAIRE(vect_residu,vect_residu,n,residu)
		residu = DSQRT(residu)
		! Remise a jour de uk
		uk = u
	ENDDO
	! On incremente le nombre d'iterations
	compteur = compteur + 1
ENDDO
CALL CPU_TIME(t_fin)
WRITE(*,*) 'Le nombre d iterations pour la methode d Euler est :'
WRITE(*,*) compteur
WRITE(*,*) 'La valeur du residu est :'
WRITE(*,*) residu
WRITE(*,*) 'Temps execution :', t_fin - t_debut, 'secondes'
!
END SUBROUTINE
